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Abstract. When a stellar-mass compact object is captured by a supermassive black hole 
located in a galactic centre, the system losses energy and angular momentum by the emission of 
gravitational waves. Subsequently, the stellar compact object evolves inspiraling until plunging 
onto the massive black hole. These EMRI systems are expected to be one of the main sources 
of gravitational waves for the future space-based Laser Interferometer Space Antenna (LISA) . 
However, the detection of EMRI signals will require of very accurate theoretical templates taking 
into account the gravitational self-force, which is the responsible of the stellar-compact object 
inspiral. Due to its potential applicability on EMRIs, the obtention of an efficient method to 
compute the scalar self-force acting on a point-like particle orbiting around a massive black hole 
is being object of increasing interest. We present here a review of our time-domain numerical 
technique to compute the self-force acting on a point-like particle and we show its suitability to 
deal with both circular and eccentric orbits. 



1. Introduction 

It is known that the majority of the galactic nuclei harbour a (super)massive black hole (MBH). 
Such a MBH can gravitationally capture a stellar-mass compact object (SCO), with mass in 
the range m = 1 — 50M Q , which eventually will perform a slow inspiral due to the emission 
of gravitational radiation until it plunges onto the MBH. These events are usually known as 
Extreme-Mass-Ratio Inspirals. EMRIs with MBHs with masses in the range M = 10 4 — 10 7 M Q 
are on the main targets of the future space-based gravitational-wave observatory LISA. One main 
issue is that EMRI signals will be buried in the LISA instrumental noise and the gravitational- 
wave foreground produced by compact galactic binaries inside the LISA band. In order to pull 
out these signals from the noise, we will have to cross-correlate the LISA data stream with a 
bank of EMRI waveform templates. These templates should be accurate enough to capture the 
physical relevant parameters of the system. For this reason, we need to include very precisely 
the gravitational back-reaction of the SCO on its own trajectory when modelling EMRIs. Due 
to the extreme mass ratios involved, an EMRI can be modeled within the framework of black 
hole perturbation theory, where the SCO is pictured as a point particle moving in the spacetime 
geometry of the MBH and the gravitational back reaction is described as the action of a local 
force, the self-force. 

As a testbed for the development of numerical techniques to compute the self-force, we use a 
simplified EMRI system, in which the SCO is a scalar charged particle, with charge q, which orbits 
a non-rotating Schwarzschild MBH [lj and whose inspiral is driven by a self-force generated by 
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Figure 1. One-dimensional spatial grid for circular orbits (left) and eccentric orbits (right). 



a scalar field, <£. This simplified model contains all the essential ingredients of the gravitational 
problem. The scalar field <3? generated by the particle satisfies a wave equation in the MBH 
geometry, and the scalar self-force acting on the particle is given by: = q{g^ v ' -\-u^u v ) (V v $)| 7 , 
where = dz^/dr is the particle's velocity , and r is proper time. 

The spherical symmetry of the MBH spacetime implies that the retarded field can be 
decomposed into decoupled scalar spherical harmonics modes: 

oo i 

*(*) = E E & m {t,r)Y im {QA) ■ (1) 

1=0 m=-l 

However, due to the singular character of the point particle, the sum over all the retarded 
field modes diverges at the particle location and, in order to obtain a self-force with a truly 
physical meaning, the field modes must be regularised. In order to do so, we use the mode 
sum regularisation scheme [21 El El El El E] , which gives an analytical expression for the singular 
contribution of the retarded field, $ s , at the particle location. Then, by subtracting the singular 
field modes from the retarded ones we obtain a smooth and differentiable field, $ R = — <3? s , 
which generates the scalar self-force = q{g^ a + u IM u a )'Va& R - Here, we review our work on 
the development of new techniques for the numerical computation of the retarded field, We 
also show the results we have obtained for the self-force after applying the mode-sum scheme. 



2. Avoiding singularities: The Particle-without-Particle Formulation 

In order to avoid the presence of the singularity associated with the point particle and all 
the related issues related with its numerical resolution we have proposed a new time-domain 
computational technique: the Particle-without-Particle (PwP) scheme. This technique was 
introduced in |8l [9] for the case of a charged scalar point particle orbiting in a circular orbit 
around a non-rotating MBH. In a subsequent work [10] . the procedure was adapted to the case 
of eccentric orbits (see also [11]) . 

The main feature of our numerical scheme consists in splitting the one-dimensional spatial 
computational domain, parameterized in terms of the radial tortoise coordinate r* , into two 
regions or subdomains, one at the left and one at the right of the particle: r* G [—00,00] = 
[— 00, r*] U [r*,oo], where r* denotes the particle radial location. Based on this division of 
the computational domain, all the variables associated with our problem can be written as: 
U(t,r*) = U_(t,r*)Q(r*(t) — r*) + U + (t,r*)Q(r*—r*(t)) , where is the Heaviside step function. 
In this way, by defining the jump of the variables as [U] = lim U + (t,r*) — lim U_(t,r*) 




r* [M] r* [M] r * [M] 



Figure 2. Snapshots of the evolution of the scalar field variable = r <I> (left), its time derivative 
dt^f (center), and its radial derivative d r ^ (right), for a particle in an eccentric orbit with orbital 
parameters (e = 0.2, p = 7 M). The internal frames show the ability of our method to resolve the 
different variables at the particle location, despite of the jumps in the time and radial derivatives. 
The snapshots have been taken after an evolution time i fina i — initial = 860M. 



(see |12[ [{JJ [TO]) we can obtain the jumps (matching conditions) at the particle location. Using 
this information we evolve the 1 + 1 wave-type equations for each mode inside each subdomain 
(which are homogeneous equations thanks to the PwP scheme). The solutions obtained are 
communicated through the boundaries by enforcing the analytic jumps. We have developed 
two practical methods for enforcing the jumps: (i) The penalty method, and (ii) the direct 
communication of the characteristic fields (see |10|). This is the essence of the PwP scheme, i.e. 
replacing the particle by a set of boundary conditions. 

In our implementation of the PwP scheme, the spatial domain is discretized using the 
PseudoSpectral Collocation (PSC) method. In this way, each subdomain is discretized 
independently using a Lobatto-Chebyshev grid and the variables are expanded in a basis of 
Chebyshev polynomials (see |13| [9]). One of the most interesting properties of the PSC 
method is that it provides exponential convergence for smooth functions. The implementation 
for the circular case has no special difficulties as the particle is always at the same radial 
position (see Figure [l) and the domain structure is static and given by $7 = |J^ =1 O a , where 

= [ r a,L; r a,ii] j an d the particle is placed at the interface of two subdomains, say Q a and r2 a+1 , 
with r* = r* K = r* +1 L = const. However in the eccentric case, the radial position of the particle 
changes with time. In this case, we keep the particle at the interface by using a time-dependent 
mapping between the spectral domain, [—1,1], and the corresponding physical subdomains, so 
that the particle is always mapped to the boundaries of the spectral domain (see Figure [T]) . 



3. Discussion 

Here we have summarized our work on the development of new time-domain scheme for the 
computation of the self-force [TO], the PwP scheme. We have shown it is a suitable scheme by 
presenting results on the simplified case of a scalar point particle in circular and eccentric orbits 
around a non-rotating MBH. In Figure [2] we show snapshots of the evolution of the different 
variables used in our implementation both for circular and eccentric orbits. This figure shows 
that the computational scheme can resolve in a very smooth way the retarded field and its 
derivatives at the particle location. In particular it can resolve jumps in the derivatives of the 
scalar field, due to the particle (this is how the particle appears in our PwP scheme). Then, by 
applying the mode-sum regularization scheme we obtain the self-force components. In Table [T] we 
show the self-force results obtained for the particle in a circular orbit, (e = 0.0, p = 6), and also 



Table 1. Self- Force components obtained for l max = 20 and 24 subdomains 



Circular Orbit Case (e = 0.0, p = 6 M) 



d t $ R 




d r $ R 


d^ R 




3.608899 


10- 4 


1.677117- 10' 4 


-5.304114 


10~ 3 








Eccentric Orbit Case 


(e = 0.2, p 


= 7 M) 


d t $ R 




d r <5> R 


d^ R 




4.278105 


io- 4 


1.946249 • 10" 4 


-5.961118 


10~ 3 



in a eccentric orbit with orbital parameters: (e = 0.2, p = 7M). In the circular case, the radial 
component of the self-force, the only one that needs to be regularized, agrees within a 1 • 10 % 
with the results obtained by [14] in the frequency domain. In Table [T] we also present results 
obtained for the eccentric orbit just mentioned. These results have been obtained by employing 
between 20 to 24 subdomains and 50 collocation points per domain. The average time for a full 
self-force calculation with £ m ax = 20 (which involves the calculation of 231 harmonic modes) in 
a computer with two Quad-Core Intel Xeon processors at 2.27GHz is around 10 to 15 minutes. 
See |11] for a detailed description of how to tune the method for performance and efficiency In 
another recent work [15J we discussed the question of initial conditions and how to avoid the 
presence of spurious solutions that can affect self-force computations. 
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